Dynamical and Statistical Criticality in a Model of Neural Tissue 
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For the nervous system to work at all, a delicate balance of excitation and inhibition must be 
achieved. However, when such a balance is sought by global strategies, only few modes remain 
balanced close to instability, and all other modes are strongly stable. Here we present a simple 
model of neural tissue in which this balance is sought locally by neurons following 'anti-Hebbian' 
behavior: all degrees of freedom achieve a close balance of excitation and inhibition and become 
"critical" in the dynamical sense. At long timescales, the modes of our model oscillate around 
the instability line, so an extremely complex "breakout" dynamics ensues in which different modes 
of the system oscillate between prominence and extinction. We show the system develops various 
anomalous statistical behaviours and hence becomes self-organized critical in the statistical sense. 

PACS numbers: 87.10.+e, 05.20.-y, 89.70.+C 



Dynamical systems theory holds that systems of inter- 
est should be structurally stable: their behavior should 
not drastically change with small perturbations of the 
defining dynamics [l| . Thus high-order criticality, the si- 
multaneous presence of several critical features such as 
Hopf bifurcations, is not expected to be ever observed 
in a natural system. However natural systems lacking 
such structural stability are not infrequent: within neu- 
roscience examples include dynamically critical systems 
such as line attractors [2] in motor control 0] and deci- 
sion making and self-tuned Hopf bifurcations in the 
auditory periphery [H| and olfactory system Q. There 
are also examples in neuroscience of statistical critical- 
ity : spontaneous heavy-tailed or scale-free fluctua- 
tions typical of critical phase transitions, such as neu- 
ronal avalanches in cortical slices anomalous correla- 
tions in the retina Q and in functional imaging [lfj| , and 
models based on simulations of the highly non-linear dy- 
namics of spiking elements, display avalanche-like statis- 
tical criticality 11 , 13| ■ There is no real understanding of 
a relation between these different concepts of criticality; 
developed turbulence, a well-studied example, displays 
both statistical criticality [l]| and dynamical criticality 
(extensive number of zero Lyapunovs [3]); but a rela- 
tionship between them is far from clear. 

We present a simple model of neural tissue, an anti- 
Hebbian network which constantly forgets; this network 
spontaneously poises itself at a dynamically critical state 
in which an extensive number of degrees of freedom ap- 
proach Hopf bifurcations, becoming arbitrarily sensitive 
to external perturbations. As the dynamics controlling 
this state has itself a marginal fixed point, the eigenvalues 
do not converge but fluctuate, close to the imaginary axis; 
when they become slightly unstable, the corresponding 
mode "breaks out" and becomes more prominent, and as 
they become slightly stable the mode slowly damps out. 
This breakout dynamics displays avalanche-like activity 



bursts whose sizes are power-law distributed. Within 
these epochs the neurons of our model are slightly corre- 
lated; yet, as the number of small but significant correla- 
tions is high, the model has strongly correlated network 
states @. This system is, on the short time-scale, sensi- 
tive in bulk to any outside input, even if applied only to a 
small subset of the neurons; however, it does not learn. In 
fact, being anti-Hebbian, it constantly forgets. We show 
that we can enrich the dynamics adding, to the term 
which is anti-Hebbian respect to regular correlations, an- 
other term "positively" Hebbian to directed correlations, 
i.e., those causal in the sense of Granger [lj|. Then the 
network may learn "predictable" stimuli, yet will stay 
unable to learn noise, and will display timing-dependent 
synaptic changes reminiscent of spike-timing dependent 
plasticity (STDP, 

We now present our model. The activities of a set of 
neurons, encoded in the vector x, evolve under the synap- 
tic connectivity matrix A; meanwhile A itself evolves, at 
a slower pace a, under an anti-Hebbian rule. 
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A = a(I - xx 1 



(1) 



(2) 



where the matrix A encodes the synaptic connections, a 
is the speed of synaptic evolution, assumed slow, and I 
the identity matrix. Inputs neuronal noise and 
nonlinear limiting terms such as x 3 would normally be 
added to the RHS of eq (1) but that shall not be necessary 
for now. From eq. (2) the matrix A would stop evolving 
when the components of x have unit variance and are 
uncorrelated to one another. 

The evolution of this system is surprisingly complex 
and generates several different timescales, as shown in 
Figures 1 and 2. For a random initial A, first, the eigen- 
mode e having the eigenvalue with the largest positive 
real part starts to diverge, and as it does so, incurs a 
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FIG. 1: Relaxation of the real parts of the eigenvalues of A. 
For clarity of illustration, N — 20. At short times (s=s 1) all 
eigenvalues with positive real parts relax to having negative 
real parts; they typically overshoot and flip sign in doing so. 
On a scale given by a = 1CP 3 , all real parts relax to the 
vicinity of the real axis. Beyond this scale, all eigenvalues 
fluctuate around the real axis. 




FIG. 2: Zooming in the rightmost portion of Fig. [T] Starting 
from a fixed point of Eq. @, i.e., an antisymmetric matrix, 
the eigenvalues fluctuate around the instability line with a 
timescale m \fa 



large penalty A rs —aee T . This happens on a timescale 
of order 1 (Figure 1). After all eigenvalues are to the left 
of the imaginary axis, a second dynamical regime ensues 
in which the real part of eigenvalues increases at a rate a 
until the real part approaches zero. Finally, the eigenval- 
ues migrating to a strip around the imaginary axis, but 
instead of relaxing to this value, they oscillate around 
their equilibrium positions (Figure 2) 

It is illustrative to use a long-time approximation 
(q << 1) when we add a noise source to each neuron: 

x = Ax + f (t) (3) 

where £ is white noise, assumed small: (£i(i)£j( s )) = 
2kT5ijS(t - s) with kT < 1. First we write Eq. © 



in the basis of the eigenvectors of A, where each com- 
ponent becomes an Ornstein-Uhlenbeck process (OUP) 
with (complex) decay rate A$, and the amplitude of the 
OUP becomes -fcT/5R(A) when U(X) < 0, and divergent 
otherwise. Then we compute (xx T |A), the expectation 
value of the correlation of the x under the assumption 
that A is constant, and use this value in Eq. 



A = a(I- 



(4) 



Diagonalize AV = VA where V are the right eigenvectors 
and A the diagonal matrix of eigenvalues Aj. Define 

A = a(I - VBV) 

The elements of B are given by 



Bij - 



A* + A* 



whose diagonal elements are 1/23JA; VBV is related to 
the inverse of the matrix VtftAV^ 1 (which we may call 
the "real part" of A). If A had orthogonal eigenvectors, 
then V^ 1 — V and hence B would be the diagonal matrix 
having 1/(23?A) in the diagonal, from where in the steady 
state A = we would obtain kT = 9?A. 

Because A is symmetric, A evolves as A(0)+S(i) where 
S is a symmetric matrix; the antisymmetric part of A 
is a constant of the motion. The evolution annihilates 
the symmetric part of A until the only piece left is the 
identity matrix multiplied by —kT, which is needed so the 
amplitudes of the OUPs isl. Therefore the fixed points of 
Eq. ©has the form A* = ±(A(0)- A T (0))-fcTI, i.e., the 
antisymmetric component of the starting matrix, minus 
kT times the identity. 

However, this fixed point turns out to be only 
marginally stable. This is easily verified: take a time 
derivative of Eq. [2] 



A 



dt 



-a — xx 

dt 



-a(xx T + xx T ) (5) 



Using x = Ax + and then reasserting that at the 
fixed point I = xx T 



A= -a(A + A T ) 



(6) 



from where we see the precise kind of marginal stability 
in question: the symetric component of the matrix A fol- 
lows an undamped harmonic oscillator equation. So this 
fixed point is itself a multidimensional Hopf bifurcation. 
The real parts of the eigenvalues of A oscillate around 
the instability boundary with a frequency y/2a. There is 
therefore a new timescale given by (Fig. [2]) . White 

noise such as in Eq. ([3]) may decohere the dynamics a 
bit, but is not necessary for driving it; even if the noise is 
turned off entirely, the system continues to oscillate. The 
oscillation frequency y2a is approximately the geomet- 
ric mean between the neuronal oscillation timescale (in 
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FIG. 3: Statistical criticality in our model. Top row, globally coupled (unrestricted A). Bottom row, nodes arranged in one 
dimension with periodic boundary conditions; only entries of A up to third nearest neighbour are allowed to be nonzero. First 
column, a display of the spationtemporal dynamics. Second column: the distribution of the number of simultaneously active 
units in the dynamics (blue) and in surrogate data (red); compare to [9(. Third column, sizes of avalanches (blue), vs. surrogate 
data (red); note in the ID case the power-law distribution of avalanche sizes, while the globally coupled (oo-D) case shows a 
piece of a power-law followed by a large lump of rather large avalanches (as clearly visible in the spatiotemporal plot). Fourth 
column, marginal distribution of the values of x (invariant under surrogation) . 



this Letter, w 1) and the synaptic update timescale a. 
In a real neuron, the oscillatory timescale is bound to be 
in the 10-120 Hz frequency bands, while the synaptic up- 
date timescale would be in the several minutes. The ge- 
ometric mean between these, marking the scale in which 
any given mode would spontaneously activate and deac- 
tivate, would be in the seconds range, bridging these two 
scales; this timescale marks the ability of our model to 
"switch task" and corresponds roughly to the timescale 
of thought. It may seem counter- intuitive that such a 
fast timescale would be a consequence of slower synap- 
tic update, unless one realizes that only in conjunction 
with the population dynamics close to a dynamically crit- 
ical state, minute alterations in synaptic strength across 
a population of neurons may change dynamical behavior 
in a mere fraction of the time required to swing a single 
synapse from low to high strength. 

The largest oscillations in Fig. [5] show an asymmetry 
that is interesting to explain. Consider our equations in 
the simple noiseless case of a single degree of freedom: 



x = Xx 

A = u(l-x 2 ) 



(7) 
(8) 



The system has two marginally stable fixed points [x = 
±1), responsible for the oscillations. Excursions on the 
first quadrant, left of a; = 1, produce runaway behav- 
ior of x and increase in A, until the quadratic term in 



Eq. ([8]) kicks in before A turns negative, leading to dif- 
ferent slopes in the up- and down-swings. This one- 
dimensional example also illustrates the origin of the sta- 
tistical criticality shown in Fig. [3] in the absence of driv- 
ing noise: the marginal stability of the modes and the 
mismatch of the dynamical and learning time-scales in- 
teract to create a long-tailed distribution of activity and 
avalanche sizes. 

In an attractor neural net, such as a Hopfield net, the 
antisymmetric components of A arc cither null or small, 
and learning is carried out by using a Hebbian rule, which 
then encodes the learned objects in the symmetric part 
of A. In our case, our anti-Hebbian dynamics takes con- 
trol of the symmetric part of A and uses it to create the 
critical, highly plastic state we have described. The anti- 
symmetric component of A is evidently untouched by Eq. 
(2), an invariant of the motion, and is the only degree of 
freedom available for learning in our system. The xx T 
term is an instantaneous density of correlation, which eq 
(2) integrates in time due to the smallness of a; let us 
write it as 



(xx T ), 



S(s)xi(t + s)xj(t — s)ds 



(9) 



An antisymmetric analog of this correlation density in 
the RHS of Eq. (2) would be given by partial directed 
correlations, correlation functions which attempt to iso- 
late influences between time-series embodying Granger 
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causality [l5|. Such correlation functions are obtained 
through a kernel which is antisymmetric in time, causing 
the correlation density to become antisymmetric in the 
neuron indexes; the simplest analog of Eq. [5]would be to 
replace the S(s) by the Hilbert transform 1/s to obtain 



C, 



Xi(t + s)xj(t — s) 



ds 



(10) 



The Hilbert transform kernel 1/s is divergent at both 
short and long timescales, and should be both ultraviolet 
and infrared cutoff according to the fastest and longest 
timescales which the system can use for its evaluation; 
the fast timescale controls the transition between increas- 
ing synaptic strength when the presynaptic neuron leads 
the postsynaptic one, to decreasing it in the opposite 
case, and reflects the accuracy with which the system 
can compute simultaneity. The slow timescale controls 
how much memory is kept of previous activity, i.e., over 
which range time intervals the pre- and post- synaptic 
activities are evaluated. When the Hilbert kernel is cut- 
off according to these two timescales, the synaptic rule 
left looks precisely like STDP 

Our results suggest a different light in which multi- 
electrode recordings may be fruitfully looked at. Our 
model proposes a view of neuronal tissue as showing coex- 
istence and superposition of different modes of neuronal 
activity, which can be simultaneously long-lived in terms 
of the timescales of electrical activity, yet extremely fast 
in terms of synaptic update timescales. The fundamental 
distinguishing factor between each of the activated modes 
is the different phase relationship of each neuron with re- 
spect to the underlying oscillation. Analysis methods 
which aim to tease apart these epochs of behavior can be 
devised by understanding that the different modes are 
distinguished from one another by looking at the activ- 
ity of single units in a coordinate frame constructed from 
the activity of the other units, rather than of external 
references. Finally it is worthwhile to remark that, since 
the dynamics consists of the activation and deactivation 
of modes of behaviour given by eigenvectors which are in 
general delocalized, the dynamics of our net is resilient to 
stochasticity or even failure in individual units. Detailed 
analysis of this resilience shall be carried out elsewhere. 
Similarly, because the dynamical modes are delocalized, 
the system is sensitive to the topological structure of the 
underlying network on scales much longer than individual 
connections or plaquettes. This extended spatial sensi- 
tivity mirrors the extended temporal behaviour discussed 
above and will be explored elsewhere. 

We have presented a simple model of "neural tissue" , in 
which an underlying anti-Hebbian dynamics permits the 
system to use the symmetric components of its synaptic 
connectivity to poise itself at a dynamically critical state 
and becomes infinitely susceptible to input which, once 
applied, can reverberate for long times. In the absence 
of inputs, this state evolves by the eigenvalues oscillating 



around the stability line, so different modes (eigenvec- 
tors) break out and then extinguish haphazardly, with 
a timescale which bridges the electrical and synaptic 
timescales. We have shown that learning can be encoded 
in the anti-symmetric component of the synaptic con- 
nectivity, driven by a term anti-symmetric both in space 
as well as time — only inputs which are Granger-causal 
and time-symmetry broken can be learned by this sys- 
tem. We have analyzed the statistics of our system to 
show that it can generate anomalous, heavy tailed dis- 
tributions, as well as power-law avalanches, showing ex- 
plicitly a connection between criticality in the dynamical 
and statistical senses. Finally, our model is intended to 
provide a scaffold to explore the implications of the re- 
verberating circuit theory introduced by Lorente de No 
and furthered by Lashley and Hebb [l7| which, for all 
their influence in physiology and behavior science, have 
not found consistent formal expressions. Supported in 
part by MCI project CGL2008-06245-C02-02 and CSIC 
intramural project HIELOCRIS (OP). 
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APPENDIX (SUPPLEMENTARY MATERIALS) 

Following the definitions of Eq. (|3|) , let us decompose 
A as AV = VA where V are the right eigenvectors and A 
the diagonal matrix of eigenvalues A.;. Let us note that 
we can define a real matrix all of whose eigenvalues are 
purely imaginary by doing H = A — V^kh^V^ 1 , this will 
be useful later. 

Multiplying (3) by V^ 1 on the left, and defining y(t) — 
V-^it) and r)(t) = 



y = Ay + t]{t) 



which has the solution 



y(t) = f 



(11) 



(12) 



from where 



(XX) 



summed over jklmn. As the integral becomes 



,(Ak+A^)(t-i') 



At, + A* 



we eventually get 



further simplified to 



\ - ^3 V jk V lk* V *ml 

4^ A,- + A 



jki 



or in other words, 



Dj 



(v-'v-ih 

Aj + A* 



(16) 



(17) 



(18) 



(19) 



Ve K ^-^V- l 2kTI5{t' - t")V e^-^Vdt'dt" 

(13) 

thus kT/2 times 

Ve^-^V-W^e^-'^Vdt'dt" (14) 

In coordinates 

/ v ij e x ^- t 'H jk v^v^*e x ^- t '^ mn v; n dt l (15) 



and notice that it becomes 1/25JA on the diagonal, 



A = ot(I - VBV) 



(20) 



where VBV is some kind of weird inverse of "Real(A)", 
namely, VSRAV -1 ; also we note that if A had orthogonal 
eigenvectors, then V^ 1 — V and hence B would be the 
diagonal matrix having 1/(23?A) in the diagonal, from 
where we'd get that for the steady state A = we would 
get kT = SKA. 



